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INTRODUCTION 


The  flow  and  heat  transfer  in  the  projectile  launching  tube  of  a 
weapon  is  typically  a  complicated  two-phase  flow  where  combustion 
products  are  mixed  with  unburned  propellant  grains.  A  detailed  calcu¬ 
lation  of  the  flow  field  in  the  gun  tube  would  provide  important 
information  such  as  local  transient  heat  transfer  rates  and  propellant 
burning  characteristics.  This  information  would  contribute  to  the 
understanding  and  solution  of  problems  associated  with  gun  barrel 
erosion  and  catastrophic  gun  failures. 

The  most  sophisticated  modeling  of  flow  phenomena  in  guns  prior  to 
the  present  work  has  been  limited  to  quasi-one-dimensional  inviscid 
two-phase  flow  analyses  of  the  propellant  combustion  process  (e.g.. 
Refs.  1-6)  and  to  time-dependent  boundary  layer  analyses  applied  to  the 


1.  Gough,  P.S.:  Numerical  Analysis  of  a  Two-Phase  Flow  with  Explicit 
Internal  Boundaries.  IHCR  77-5,  Naval  Ordnance  Station,  Indian 
Head,  MD,  April  1977. 

2.  Koo,  J.H.  and  Kuo,  K.K. :  Transient  Combustion  in  Granular  Propel¬ 
lant  Beds.  Part  1:  Theoretical  Modeling  and  Numerical  Solution  of 
Transient  Combustion  Processes  in  Mobile  Granular  Propellant  Beds. 
BRL  CR.-346,  U.S.  Army  Ballistic  Research  Laboratory,  Aberdeen 
Proving  Ground,  MD,  August  1977.  (AD  #A044998) 

3.  Kuo,  K.K.,  Koo,  J.H.,  Davis,  T.R.  and  Coates,  G.R.:  Transient 
Combustion  in  Mobile  Gas-Permeable  Propellants.  Acta  Astronautica , 
Vol .  3,  1976,  pp.  573-591. 

4.  Fisher,  E.B.,  Graves,  K.W. ,  and  Trippe,  A.P.:  Application  of  a 
Flame  Spread  Model  to  Design  Problems  in  the  155  mm  Propelling 
Charge.  12th  JANNAF  Combustion  Meeting,  CPIA  Publication  273, 

Vol.  I,  December  1975,  p.  199. 

5.  Krier,  H.  ,  Raian,  S.,  and  VanTassell,  W.  :  Flame  Spreading  and 
Combustion  in  Packed  Beds  of  Propellant  Grains.  AIAA  Journal, 

Vol.  14,  No.  3,  March  1976,  p.  301. 

6.  Krier,  H.  and  Gokhale,  S.S.:  Modeling  of  Convective  Mode  Combustion 
Through  Granulated  Propellant  to  Predict  Detonation  Transition. 

AIAA  J.,  Vol.  16,  No.  2,  1978,  pp.  177-183. 
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flow  of  propellant  gases  in  a  gun  barrel  (Refs.  7  and  8).  The  boundary 
layer  procedures  suffer  from  the  shortcoming  that  the  starting  condi¬ 
tions  near  the  projectile  base  are  not  well  defined,  and  according  to 
conventional  boundary  layer  theory  the  heat  flux  near  the  base 
approaches  infinity  because  the  base  is  a  singular  point  (Ref.  8). 
Furthermore,  the  validity  of  the  boundary  layer  approximations  is 
questionable  at  both  the  breech  end  and  the  projectile  base  region,  and 
even  the  most  sophisticated  boundary  layer  analysis  presently  used  for 
gun  barrel  problems,  e.g..  Refs.  7  and  8,  did  not  consider  the  two-phase 
flow  aspects  of  the  propellant  combustion  process.  The  significant 
features  of  the  two-phase  flow  interior  ballistics  codes  (Refs.  1-6) 
were  reviewed  recently  by  Kuo  (Ref.  9).  The  main  objection  to  these 
analyses  (Refs.  1-6)  would  seem  to  be  the  presumption  of  quasi-one- 
dimensional  flow  and  the  attempt  to  predict  heat  transfer  to  the  barrel 
using  rather  simple  unsteady  boundary  layer  models  or  correlation 
formulas . 

Under  the  present  effort  a  mathematical  model  of  a  two-phase,  two- 
dimensional  flow  was  developed  and  a  computer  code  has  been  constructed 
for  the  numerical  solution  of  the  equations  resulting  from  this  mathe¬ 
matical  model.  The  model  developed  consists  of  the  governing  equations 
for  an  axisymmetric ,  two-phase  flow  in  a  gun  tube  with  a  rotating 
projectile,  and  a  system  of  constitutive  relations  describing  the 
molecular  viscosity  and  thermal  conductivity,  turbulence  length  scale, 


7 .  Anderson,  L.W. ,  Bartlett ,  E.P . ,  Dahm,  T. J .  and  Kendall ,  R.M. : 
Numerical  Solution  of  the  Nonsteady  Boundary  Layer  Equations  with 
Application  to  Convective  Heat  Transfer  in  Guns.  Aerotherm  Report 
No.  70-22,  Aerotherm  Corp . ,  October  1970. 

8.  Bartlett,  E.P.,  Anderson,  L.W.,  and  Kendall,  R.M. :  Time-Dependent 
Boundary  Layers  with  Application  to  Gun  Barrel  Heat  Transfer. 
Proceedings  12th  Heat  Transfer  and  Fluid  Mechanics  Institute, 
Stanford  Univ.  Press,  1972,  p.  262. 

9.  Kuo,  K.K. :  A  Summary  of  the  JANNAF  Workshop  on  "Theoretical 
Modeling  and  Experimental  Measurements  of  the  Combustion  and  Fluid 
Flow  Processes  in  Gun  Propellant  Charges".  13th  JANNAF  Combustion 
Meeting,  CPIA  Publication  281,  Vol .  I,  December  1976,  p.  213. 
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gas  equation  of  state,  intergranular  stress,  interphase  drag,  interphase 
heat  transfer,  and  solid  phase  combustion.  The  governing  equations  and 
corresponding  initial  and  boundary  conditions  describe  the  firing  cycle 
beginning  with  a  fluidized  and  ignited  solid  phase,  and  ending  with  the 
projectile  exiting  the  gun  tube.  Chemical  reactions  within  the  gas 
phase  were  excluded  from  the  formulation.  An  axisymmetric  time-dependent 
adaptive  coordinate  system  for  interior  ballistics  flow  field  calcu¬ 
lations  was  developed,  and  the  projectile  and  distinct  filler  elements 
were  treated  using  a  quasi-one-dimensional  lumped  parameter  analysis. 

The  complex  nature  of  the  flow  in  the  projectile  base  region  and 
in  the  breech  end  of  the  barrel  does  not  permit  simplifying  approxi¬ 
mations  to  be  made  in  the  governing  fluid  flow  equations,  and  therefore 
in  principle  the  solution  of  the  full  Navier-Stokes  equations  is 
required,  rather  than  some  simpler  approximate  set  of  equations. 
Fortunately,  recent  developments  in  computational  fluid  dynamics  have 
made  possible  the  prediction  of  the  detailed  flow  field  in  configur¬ 
ations  such  as  a  gun  barrel  using  the  full  Navier-Stokes  equations. 

TKg  equations  and  coordinate  system  developed  under  this  effort  have 
been  incorporated  into  an  existing  three-dimensional  time-dependent 
compressible  Navier-Stokes  calculation  procedure  (the  MINT  code)  which 
was  originally  developed  under  United  States  Navy  and  Air  Force  sponsor¬ 
ship  for  other  purposes  by  staff  members  of  Scientific  Research 
Associates,  Inc.  (Refs.  10-13).  The  MINT  procedure  solves  the  governing 


10.  Briley,  W.R. ,  and  McDonald,  H. :  An  Implicit  Numerical  Method  for 
the  Multidimensional  Compressible  Navier-Stokes  Equations.  United 
Aircraft  Research  Laboratories  Report  M911363-6,  November  1973. 

11.  Briley,  W.R.,  McDonald,  H. ,  and  Gibeling,  H.J.:  Solution  of  the 
Multidimensional  Compressible  Navier-Stokes  Equations  by  a 
Generalized  Implicit  Method.  United  Technologies  Research  Center 
Report  R75-911363-15 ,  January  1976. 

12.  Briley,  W.R.,  and  McDonald,  H. :  Solution  of  the  Multidimensional 
Compressible  Navier-Stokes  Equations  by  a  Generalized  Implicit 
Method.  J.Comp.  Physics,  Vol.  24,  No.  4,  1977,  p.  372. 

13.  Gibeling,  H.J.,  McDonald,  H.,  and  Briley,  W.R.:  Development  of  a 
Three-Dimensional  Combustor  Flow  Analysis.  AFAPL-TR-75-59 ,  Vol.  I, 
July  1975  and  Volume  II,  October  1976. 
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equations  using  a  consistently-split,  linearized,  block-implicit 
numerical  scheme  (Ref.  14).  The  resulting  computer  code  will  be 
designated  as  the  MINT-G  code  herein. 


14.  Briley,  W.  R.  and  McDonald,  H. :  On  the  Structure  and  Use  of  Linear¬ 
ized  Block  ADI  and  Related  Schemes.  SRA  Report  R78-3A,  to  appear 
in  J.Comp.  Physics,  1979. 
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THEORETICAL  ANALYSIS 


Approach 

The  governing  equations  for  a  two-phase  two-dimensional  flow  in  a 
gun  tube  are  presented  below.  The  provision  for  a  rotating  projectile 
is  considered  by  solving  the  azimuthal  momentum  conservation  equation 
with  the  appropriate  boundary  conditions  at  the  projectile  base.  The 
governing  equations  may  be  obtained  by  employing  either  the  time¬ 
averaging  procedure  utilized  by  Ishii  (Ref.  15)  or  the  formal  averaging 
approach  used  by  Gough  (e.g.,  Refs.  16,  17)  or  Gough  and  Zwarts  (Ref.  18) 
In  the  present  derivation,  the  averaging  procedure  of  Gough  (Ref.  16) 
has  been  selected  because  of  its  notational  convenience;  however, 
extensive  reference  to  the  work  of  Ishii  (Ref.  15)  has  been  made  in 
order  to  verify  the  results  obtained.  In  the  following  analysis,  a 
gas-solid  mixture  is  assumed  with  a  constant  solid  phase  density,  p  . 
Numerous  assumptions  and  approximations  are  required  in  order  to  formu¬ 
late  a  tractable  problem.  Most  of  the  required  assumptions  have  been 
stated  previously  by  Gough  (e.g..  Ref.  1),  and  those  necessary  in  the 
present  work  are: 

(1)  The  gas  and  solid  phases  occupy  separate  complementary  regions,  and 
within  each  region  the  material  may  be  treated  as  a  homogeneous  continuum 

(2)  The  flow  of  the  heterogeneous  mixture,  composed  of  the  two  inter¬ 
acting  continua,  can  be  described  by  appropriately  defined  averages  of 
the  flow  properties. 


15.  Ishii,  M. :  Thermo-Fluid  Dynamic  Theory  of  Two-Phase  Flow. 
Eyrolles,  Paris,  1975. 

16.  Gough,  P.S.:  Derivation  of  Balance  Equations  for  Heterogeneous 
Two-Phase  Flow  by  Formal  Averaging.  ARO  Workshop  on  Multiphase 
Flows,  Ballistic  Research  Laboratory,  February  1978,  pp.  71-80. 

17.  Gough,  P.S.:  The  Flow  of  a  Compressible  Gas  Through  an  Aggregate 
of  Mobile,  Reacting  Particles.  Ph.D.  Thesis,  Department  of 
Mechanical  Engineering,  McGill  University,  Montreal,  1974. 

18.  Gough,  P.S.  and  Zwarts,  F.J.:  Some  Fundamental  Aspects  of  the 
Digital  Simulation  of  Convective  Burning  in  Porous  Beds.  AIAA 
Paper  77-855,  July  1977. 
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(3)  If  solid  phase  combustion  occurs,  the  energy  deposition  is  taken 
to  be  in  the  gas  only. 

(4)  The  solid  phase  is  deformable  and  incompressible.  However,  locally 
no  relative  motion  between  the  solid  particles  is  considered.  Thus  the 
average  stress  in  the  solid  phase  is  an  isotropic  normal  stress. 

(3)  The  influence  of  solid  phase  deformation  on  the  particle  surface 
area  is  neglected,  and  the  interfacial  average  of  the  particle  velocity 
is  equal  to  the  volume  average  in  the  absence  of  burning. 

(6)  The  interphase  drag  is  determined  from  steady  state  correlations; 
the  unsteady  virtual  mass  effect  is  not  considered. 

(7)  The  interphase  heat  transfer  is  determined  from  steady  state 
correlations . 

(8)  The  Noble-Abel  equation  of  state  will  be  employed.  The  specific 
heats  (c  and  c^)  are  taken  to  be  independent  of  temperature. 

(9)  The  regression  rate  of  the  surface  of  the  burning  propellant  is  a 
function  of  the  average  gas  properties  and  the  propellant  surface  tem¬ 
perature  . 

(10)  Heat  transfer  to  the  solid  phase  is  treated  as  a  one-dimensional 
process  in  order  to  determine  the  propellant  surface  temperature. 

(11)  The  pressure  drop  at  the  gas-solid  interface  is  negligible. 

Governing  Equations 

Both  Ishii  (Ref.  15)  and  Gough  (Refs.  16,  17)  have  presented  the 

relations  for  the  average  of  time  and  space  derivatives  in  a  two-phase 

mixture.  Using  the  above  assumptions  a  system  of  partial  differential 

equations  is  obtained  containing  interface-averaged  source  terms  arising 

from  averaging  the  basic  conservation  equations  for  the  two-phase 

mixture.  A  basic  quantity  used  to  describe  a  two-phase  mixture  is  the 

porosity,  a,  i.e.,  the  ratio  of  volume  occupied  by  the  gas  phase  to  the 

total  volume.  Ishii  (Ref.  15)  introduces  several  averages  which  are 

required  in  the  present  analysis.  Gough  (Ref.  16)  introduces  a  general 

^  — * 

weighting  function  g(y-x,i-t)  which  reflects  the  influence  of  remote 
points  (y,r)  on  the  average  value  at  (x,t).  By  definition,  the  Gough 
average  gives 
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(1) 


J  Q(x,t)  dx  dt  =1 
oil  v,r 

The  porosity  is  defined  by 

a(x  -  J  g(y -x  ,  T-t)  dy  dr  (2) 

Vgas 

The  .eighting  function,  g,  plays  a  role  similar  to  the  state  density 
functions  (M^,  ,  M  )  introduced  by  Ishii  (Ref.  15,  p.  65).  The  basic 

time  average  introduced  by  Ishii  (Ref.  15,  p.  68)  is  denoted  by  a 
single  overbar  (4») ,  and  this  is  equivalent  to  Gough's  (Ref.  2) 
unnormalized  average.  The  phase  average  denoted  by  a  double  overbar 
(ijj )  is  related  to  ip  by 

=  ^  I  r  _  w 

'P  ■  —  *  —  J  gly-x,  T-t)\f/(y,T)dydr  (3) 

a  V 
vgos 

Eq.  (3)  defines  the  average  of  a  gas  property,  i p,  since  the  integral  is 
carried  out  over  the  region  occupied  by  the  gas  phase,  .  In  Ishii* s 

approach  the  equivalent  average  is  obtained  by  integrating  only  over  the 
time  interval  for  which  the  gas  phase  is  present  at  the  space  point  x. 
Finally,  the  mass  weighted  average  for  a  property  of  the  k^-phase  ^ 
is  defined  by 

r  .  .  ML  (4) 

k  K  Py 

This  average  is  also  known  as  the  Favrd  average,  hence  the  superscript 

F  is  used.  This  is  a  very  convenient  average  to  use  in  turbulent  flow 

since  density  fluctuations  may  be  eliminated  formally.  It  should  be 

_  t 

noted  that  the  quantity  is  the  partial  density  of  k  -phase  while 
p  is  the  actual  density,  so  that  the  mixture  density  is  given  by 
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(5) 


2 


2 

Z  a 

k=| 


krk 


where  =  a,  and  -  1-a. 

In  the  following  equations,  the  Favre  average  is  introduced  where 
it  is  appropriate,  and  phase  average  values  are  used  otherwise.  The 
Favr£— averaged  velocity  vector  is  written  as 

UF  =  U  (6) 


and  on  all  other  variables  (e.g.,  e,  h,  etc.)  the  superscript  F  is 
dropped  for  convenience.  The  fluctuating  component  of  any  variable  is 
denoted  with  a  superscript  prime,  All  quantities  pertaining  to  the 

solid  phase  are  denoted  by  the  subscript  p.  The  resulting  equations 
are  then 

Gas  Phase  Continuity 


d{ap ) 
Tt 


+  V  *  { ap  u ) 


(7) 


Solid  Phase  Continuity 


d(l-a) 

“5t 


+  V  •  [  ( I  -  a  )  Up] 


(8) 


where  the  mass  source,  is  due  to  propellant  burning.  Following 

Gough  (Ref .  16) , 


Tt  =  -  J  p{  u  -  u')-n  gdA  (9) 

I 

where  u  is  the  velocity  of  the  interface  between  the  phases,  I  is  the 
region  of  integration  as  defined  by  the  interphase  surface  and  time, 
and  dA  is  the  differential  element  in  I-space  (i.e.,  dA  is  an  area-time 
product).  Introducing  the  instantaneous  surface  regression  rate,  d, 
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the  interface  velocity  is 


U* 


=  u,  +  nd 


(10) 


where  n  is  the  outward  normal  from  the  gas  phase.  The  instantaneous 
interfacial  boundary  conditions  were  stated  by  Ishii  (Ref.  15, 
pp.  29-30),  and  by  Gough  (Ref.  16)  under  the  assumption  that  surface 
tension  is  zero  and  the  surface  energy  remains  constant.  These 
relations  are 


/>(u-u')-n  =  />p(up-u')'K 


(lla) 


[pu(u-u')- n}n  -  [/jpUpCUp-u1)-  np]-n  (lib) 


[ />(u  - u')( e  +  1-  u-u)  -  n-u  +  q]  •  n 

[/»p(Sp-3,)(ep  +  ^-3p-«Ip)-n  -Up  +  gp 


•  n 


(11c) 


where  q  and  q  are  heat  flux  vectors,  the  total  stress  tensor  II  is 


n  *  -pi  +  t r 


(12) 


and 


np  =  -  pn  +  r 


(13) 


Here  ]R  is  the  granular  stress  tensor  in  the  solid  phase.  The 
traditional  sign  convention  for  stress  has  been  chosen  herein,  i.e., 
negative  in  compression  and  positive  in  tension.  Using  Eqs.  (10)  and 
(lla)  in  Eq.  (9)  yields 
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(14) 


r,  =  f 


According  to  Gough  (Ref.  16)  the  interface  average  of  tp  is  defined  as 


<*>'  -  1 


J  i^gd  A  J i//gdA 

_  2 


(l-a)Sp/Vp 


(15) 


where  is  the  average  particle  surface  area  and  V  is  the  average 
particle  volume.  Hence  Eq .  (14)  becomes 


r,  -  (l-a)-^fE.  <d>' 


(16) 


where  <d>1  is  the  average  regression  rate  of  the  solid  phase.  In  Ref.  1 
Gough  has  extended  the  source  term  expression  for  F  to  include  a 
number  of  different  types  of  granular  particles,  and  if  necessary, 
this  extension  could  be  incorporated  into  the  present  formulation  at  a 
later  time. 

In  the  present  work,  Eqs.  (7-8)  would  be  solved  in  conjunction 
with  Eq-  (16)  and  a  constitutive  relation  for  <d>1. 

Gas  Phase  Momentum 


afa^U)  _  __  _ 

— ^ -  +  V  •  ( a/>  UU )  =  -  V(ap)  +  V  •  [  a  (  tt  +  t r 1  )J  +  M  (17) 


r  _  -r  1 


Solid  Phase  Momentum 


a[(i-d)/>pOp] 


at 


+  V-[(l-a)/>  U  Up] 


X7  F  N  -  rt  )  n 
*  ^  \  1  »•*  /  k 


O  .  f  l  I  - 

v  ^  \  i 


W  RD  X  I  ^ 
/  V  m  <  Up/ 
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(18) 


In  the  above  equations,  tt  and  tt  are  the  average  stress  tensor  and  the 

turbulent  stress  tensor  in  the  gas  phase,  respectively,  1R  is  the  aver- 

T 

age  granular  stress  tensor,  and  tt  is  the  solid  phase  turbulent  stress 

T  P 

tensor.  For  the  present  time  tt^  will  be  neglected  because  there  is 
insufficient  information  available  to  construct  a  constitutive  relation 
for  it.  The  gas-solid  momentum  exchange  term,  M^,  is  defined  by 


M  =  -  f  \  nit  ^  it  -  Tt*  ^  -  11 1  •  n  n  rl  A 
’I  “  J  L  r~ '  -  -  '  “ J  " 

2 


/i  nN 
\j-jj 


Assuming  the  change  in  normal  stress  at  the  interface  is 


( n p  -  IT )’  •  n  =  n  Ap1  +  (  R  -  tt)1  •  n  -  (  R  -  tt )’ •  n  (20) 


1 


and  using  Eqs.  (10)  and  (lib)  in  Eq.  (19)  ,  one  obtains 


Mj  ■=  -  J"[-(n+n')-n-(R--ir)-n-^  upd]gd 


(21) 


Further,  assume  that  n  -  it  ,  and  JR  -  U  at  the  interlace  so  that 


M, 


=  -  /  [  pK  -  n'-n  -/>pupd]  g 


dA 


(22) 


Noting  that  (Ref.  15,  p,  75) 


J n  gdA  •  - Va 


(23) 


the  expression  for  becomes 


M(*pva  +  j  (  n;*  K)  gdA  +  j  (  Up  +  ) d  go A 

Z  Z 


(24) 


Finally,  the  contribution  in  the  last  term  is  neglected,  and  the 

infprnhpQP  Hrp»a  npr  unit"  nf  cinl  nluxso  nc  Hof  i  npH  ^g  fPof  16^ 

- ^4.—  ~W - O  r - - - - *-■«-*  ^  **'-*-•  -*-W/ 


<  F  >'  ( I  -  a )  — =  -/  ( II^  n )  gdA 


(25) 


P  Z 
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so  Eq.  (24)  becomes 
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The  latter  term  is  difficult  to  model  except  for  a  dispersed  flow 
(Ref.  15,  p.  165),  hence  it  must  be  neglected  at  present.  The  turbu¬ 
lent  flow  stress  tensor  in  the  gas  phase  will  be  modeled  using  an 
isotropic  eddy  viscosity  formulation,  i.e., 

•n-T  ■=  -p  u7!7  =  2/itD|  -  -|-(  U  +  pi.  )tt  (33) 

where  k,  the  turbulence  kinetic  energy,  is  discussed  in  the  section 
on  Turbulence  Model  Equations.  The  turbulent  viscosity  must  be 
determined  using  a  suitable  turbulence  model.  The  solid  phase 
granular  stress  tensor,  H,  will  be  modeled  by  assuming  an  isotropic 
normal  stress,  i.e., 


R  £  RpI 


(34) 


hence  in  the  solid  phase  momentum  Eq .  (28) 


m  /  i  _  \ 

v  •  I  ii  -  a;  rcj 


r~»  r  / 1  _  \ « 

V 


P  J 


/  o  c  \ 


In  the  present  work,  Eqs.  (27-28)  would  be  solved  in  conjunction 

with  Eqs.  (16),  (29-31),  (33)  and  (35),  and  constitutive  relations  for 

<d>^,  <F>1  and  R  . 

P 

Gas  Phase  Energy  Equation 

In  the  present  formulation  it  is  desirable  to  write  the  energy 
equation  in  terms  of  the  mass-averaged  static  enthalpy  h  because  of 
numerical  considerations  in  solving  the  resulting  coupled  system  of 
finite  difference  equations. 


d(a^h) 

at 


+  V  •  ( a  u  h  )  =  -  V  ■  [  a  (  q  +  q  T )] 


(36) 


D  _  ^  ij .  u  * 

+  FT  (aP*  +  ®  +  aP*  +  (  E,  -  u  m(  +  q ) 
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where  <j>  is  the  mean  flow  dissipation  term  defined  in  Eq.  (A-10)  and  e 

is  turbulence  kinetic  energy  dissipation  rate.  The  mean  heat  flux  vector 
— 

q  and  the  turbulent  heat  flux  vector  q  in  a  two-phase  flow  may  be 
written  as  (Ref.  15,  p.  165) 


^  _ r  V  Cl  -  I 

q  ■  -Z[  VT-  —  (T,  -  T)]  (37) 

and 

T  r  _  Va  =  1 

qT  .  -KJ  [  VT  -  —  (T,  -  T)J  (38) 


where  <  is  the  mean  thermal  conductivity,  k  is  a  turbulent  conduct¬ 
ivity,  and  T_^  is  the  mean  temperature  at  the  interface  between  the 
phases.  For  the  present  time  will  be  taken  as  the  average  between 
the  gas  temperature  and  the  particle  surface  temperature,  i.e. 


r(T+V 


(39) 


and  T  will  be  determined  from  the  solid  phase  heat  conduction  model, 
ps 

The  effective  conductivity  will  be  modeled  using  an  effective  Prandtl 
number  obtained  from  knowledge  of  turbulent  flows  of  gases  and  gas 
mixtures ,  i.e.. 


K 


eff 


=  K  +  K 


p^eff 


Pr 


eff 


(40) 


where  the  effective  viscosity  is  the  sum  of  the  laminar  and  turbulent 
viscosities , 


Meff  ‘MVt 

A  constant  value  will  be  employed  for  the  effective  Prandtl  number 
Pr  ^  =  0.9.  Following  Gough  (Ref.  16),  it  can  be  shown  that  the 
interfacial  energy  transfer  term  in  Eq.  (36)  is 
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u  u 


E,  -0-«,  +  — 


r,  *  -p(u-up)-Va 


+  ( I  -  a )  — ^  (  U  -  Up)  •  <  F  >  +  q  •  Va 
VP 

-d-a)^  <q>j  +  r,[hcomb+-f  (U-Up)-(U-Op)] 


(42) 


where  <q>1  is  the  interfacial  average  heat  transfer  between  the  gas 
and  solid  phases,  and  h  ^  is  the  energy  released  (per  unit  mass)  due 
to  combustion  of  the  solid  propellant. 

In  the  present  work,  Eq.  (36)  would  be  solved  in  conjunction  with 

*  i  i 

Eqs.  (16),  (26)  and  (37-42),  and  constitutive  relations  for  <d>  ,  <F>  , 
and  <q>1. 

Solid  Phase  Heat  Conduction  Equation 

Since  the  solid  particle  surface  temperature  is  desired  to  deter¬ 
mine  ignition,  the  propellant  burning  rate,  and  the  rate  of  heat 
transfer  between  the  gas  and  solid  phase,  a  transient  heat  conduction 
equation  must  be  solved.  Gough  (Ref.  1)  and  Kuo,  et  al.,  (Ref.  3) 
assume  that  the  penetration  depth  of  a  thermal  wave  into  the  propellant 
grains  is  small  compared  to  the  grain  dimensions.  Then  it  is 
permissible  to  use  a  one-dimensional  approximation  (planar  for  cord, 
propellant  or  spherical  for  spherical  propellant  grains)  to  obtain  the 
particle  surface  temperature.  Following  the  motion  of  a  given  particle 
(Kuo,  et  al..  Ref.  3),  the  heat  conduction  equation  for  a  spherical 
particle  is 

/dfp\  qp  a2(r-2fp) 

\  d  t  / ~  r2  a"r2  <43> 


where  =  T^(r;  x,t)  is  the  phase-averaged  temperature  within  a 
representative  particle,  r  is  radial  coordinate  within  the  particle,  a 

is  the  thermal  diffusivity  of  the  solid  particles  [a  =  k  /p  (c  )  ] 

P  p  Hp  p/pJ* 
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and  (d/dt)~  denotes  the  Lagrangian  time  derivative  at  constant  r  within 
the  particle.  Since  the  surface  of  a  representative  burning  particle  is 
receding  in  time  it  is  desirable  to  employ  the  following  time-dependent 
transformation  for  the  particle  radial  coordinate  i: 


'■a1  osEsl 


(44) 


6Jn\  (  £  drD\dTD  aD  a  2 


v dt  yc- y-ar  V 

where  the  quantity 

drD  *  i 

p  = - =  <  d  > 

5  "  dt 


Then  Eq.  (43)  becomes, 

=  \  / 

P 


(45) 


(46) 


may  be  idenfitied  as  the  average  surface  regression  rate  for  the 
particle,  ^  0. 

The  initial  condition  for  Eq.  (45)  is 


Tp<C>t  =  0>  ‘  Tpo  (47) 


The  boundary  conditions  are 


(£'  o,t)  =  o 


at  c=o 


(48) 


*  hc(t)[f-Tps]+ qrad  +kp</,(Rs>P)  at  C-K49) 

P  ^ 

where  is  the  net  incident  radiation  heat  flux  normal  to  the 

nRAD 

particle  surface,  is  the  thermal  conductivity  of  the  solid  particles, 
and  4> ( , p )  is  the  heat  feedback  from  the  flame  identified  by  Gough 
(Ref.  1,  p.  57).  Assuming  that  the  gas  is  nearly  in  radiative  equil¬ 
ibrium  so  that  the  gas  emissivity  is  unity,  and  that  radiation  emitted 
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by  other  particles  does  not  influence  the  particle  in  question,  we 
obtain 


^rad  *  £pcr  (  “  Tps  }  (50) 

where  is  the  particle  emissivity.  Other  authors  (e.g..  Refs.  1-3) 
have  cast  Eq.  (50)  into  a  heat  transfer  coefficient  form,  so  Eq .  (49) 
becomes 

1  £ =  l,t)  *  ht(,)[f  -  v] +  kp^(Rs>P)  <51> 

where  the  total  heat  transfer  coefficient  is 


h 


t 


+  «pcr  (  T  +  Tps)  (  T2 


+  f  PS  ) 


(52) 


The  convective  heat  transfer  coefficient,  h^,  will  be  specified  via 
constitutive  relations  below.  An  expression  for  4»(Rs,p)  has  been 
presented  by  Gough  (Ref.  1)  for  a  planar  geometry  under  the  assumption 
that  the  flame  zone  surrounding  the  burning  particle  remains  quasi¬ 
steady,  and  that  the  convection  and  radiation  heat  transfer  terms 
in  Eq.  (49)  are  zero.  It  then  follows  that 


a  (TPs  TpcP 
P 

where  a  is  the  thermal  diffusivity  of  the  particles  and  T  is  the 
P  po 

undisturbed  temperature  far  from  the  particle  surface.  In  the  context 
of  spherical  particles,  would  be  taken  as  the  temperature  at  the 
center  of  the  particle.  This  procedure  should  be  sufficiently  accurate 
in  view  of  the  other  assumptions  made  in  obtaining  Eq.  (53) 

In  the  present  work,  Eqs.  (47-50,  53)  will  be  utilized.  Solution 
of  this  solid  phase  heat  conduction  model  requires  special  considera¬ 
tion  since  it  is  in  Lagrangian  form,  whereas  all  other  differential 
conservation  equations  are  in  Eulerian  form.  The  method  to  be  employed 
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in  the  present  analysis  will  be  described  in  the  section  on  Solution 
Procedure. 

Turbulence  Model  Equations 

The  introduction  of  the  turbulent  viscosity  (y  )  in  Eq .  (33) 
requires  the  use  of  a  turbulence  model  to  specify  this  quantity.  It 
was  originally  anticipated  that  a  two-equation  turbulence  transport 
model  would  be  implemented  in  conjunction  with  the  Prandt 1-Kolmogorov 
formula  for  specification  of  the  turbulent  viscosity,  i.e.. 


H- r  '  ‘I/’*"2! 


(54) 


where  k  is  the  turbulence  kinetic  energy  and  l  is  a  length  scale  of 
the  turbulence.  This  relation  follows  from  dimensional  arguments  for 
turbulent  flow  described  by  the  two  parameters,  k  and  £,  Various  forms 
of  the  two-equation  model  of  turbulence  have  been  proposed  since 

Tr_1  /  Ti  ~  £  TAN  c  ■ ..  _ _ A A  j-  1 _ t-  : TO/.  O  ^ 

f'vO J.HJO go r ov  <.rv.ei  .  ixisl  iijuuuucku  luc  i-uuLejjL  m  1744,  husl 

investigators  have  chosen  the  kinetic  energy  of  turbulence,  k,  as  their 
first  variable.  A  commonly  chosen  second  variable  has  been  the  turbu¬ 
lence  kinetic  energy  dissipation  rate,  e. 

The  appropriate  transport  equations  for  turbulence  kinetic  energy 
and  energy  dissipation  rate  valid  at  high  Reynolds  numbers  have  been 
presented  by  Launder  and  Spalding  (Ref.  20).  However ,  the  very  large 
fluid  accelerations  experienced  in  the  interior  ballistics  problem 
require  the  consideration  of  the  laminarization  of  the  turbulent  flow 
near  solid  surfaces.  There  are  two  options  available  for  modeling  the 
turbulence  near  a  wall.  In  the  first,  grid  point  resolution  normal  to 


19.  Kolmogorov,  A.N. :  Equations  of  Turbulent  Motion  of  an  Incompres¬ 
sible  Turbulent  Fluid.  IZC.  Adak.  Naut.  SSR  Ser.  Phys.  VI, 

No.  1-2,  56,  1942 . 

20.  Launder,  B.E.  and  Spalding,  D.B.:  The  Numerical  Computation  of 
Turbulent  Flows.  Computer  Methods  in  Applied  Mechanics  and 
Engineering,  Vol.  3,  1974,  p.  269. 
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the  wall  must  be  sufficient  to  define  the  viscous  sublayer,  in  which 
case  the  boundary  conditions  are  relatively  straightforward.  However, 
the  difficulty  with  this  approach  is  that  the  physics  of  low  Reynolds 
number  (transitional)  turbulence  must  be  modeled  in  a  reasonable  manner 
by  the  governing  turbulence  equations  (e.g.,  Jones  and  Launder, 

Ref.  21).  An  alternative  approach  is  to  employ  a  less  refined  mesh 
and  force  the  turbulence  variables  to  yield  values  consistent  with  a 
boundary  layer  wall  function  formulation  at  the  first  grid  point  away 
from  the  wall.  The  difficulty  with  this  approach  is  that  the  validity 
of  the  wall  function  formulation  is  questionable  under  the  rapidly 
accelerating  transient  flow  situation  present  in  the  interior  ballis¬ 
tics  problem.  Furthermore,  recent  experience  at  SRA  indicates  that 
the  wall  function  approach  will  be  inadequate  for  a  reacting  unsteady 
flow  with  moving  coordinates.  In  addition,  SRA's  experience  with  the 
k-e  turbulence  model  has  shown  it  to  be  unreliable  both  in  reacting 
flows  with  large  energy  release  and  in  complicated  transitional  flows 
where  the  viscous  sublayer  is  resolved. 

Therefore,  in  the  present  work  it  is  proposed  to  utilize  a 
turbulence  kinetic  energy  equation  in  conjunction  with  a  specified 
turbulence  length  scale  distribution.  An  equation  for  the  turbulence 
kinetic  energy  of  the  gas  phase  may  be  derived  using  the  averaging 
procedure  of  Ishii  (Ref.  15)  or  Gough  (Ref,  16).  Following  the 
derivation  of  Bradshaw  and  Ferriss  (Ref.  22),  one  obtains 


diap  k)  p. 

+  V-{ap Uk)  -  V-  (a  ~  Vk) 

u  k 


+  a  -|-(V-U)*]  -  ^-^kV-U  -  pi 


at 


(55) 


+  S, 


21.  Jones,  W.P.  and  Launder,  B.E.:  The  Prediction  of  Laminarization 
with  a  Two-Equation  Model  of  Turbulence.  Int.  J.  Heat  Mass 
Transfer,  Vol.  15,  1972,  p.  301. 

22.  Bradshaw,  P.  and  Ferriss,  D.H.  :  Calculation  of  Boundary-Layer 
Development  Using  the  Turbulent  Energy  Equation:  Compressible 
Flow  on  Adiabatic  Walls.  J.  Fluid  Mechanics,  Vol.  46,  Part  1, 
1971,  pp.  83-110. 
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where  k  is  defined  as 


k  — ♦/  —*i 

=  —  U  •  U 


(56) 


and  e  is  the  turbulence  energy  dissipation  rate.  The  interfacial 
transfer  term,  S^,  is 

Sk  =  -  /  [  d-  ^(iT'.u'Hu-u')  +  p'u']  -ngdA 


(57) 


If  the  pressure-velocity  correlation  is  neglected  and  an  average  value 
for  l/2(u*  •  u’)  is  assumed  to  be  k^  at  the  gas-solid  interface,  S, 
become  s 


ps 


Sk  “  kpsr, 


(58) 


Evidently,  this  term  represents  the  production  of  turbulence  kinetic 

energy  in  the  gas  phase  due  to  gasification  of  the  solid  particles. 

However,  it  is  not  known  how  to  specify  k  at  the  present  time. 

ps 

Using  dimensional  arguments  the  Prandtl-Kolmogorov  formula, 

Eq .  (54)  may  be  written  as 


H-  T  =  c 


H-  € 


(59) 


and  the  dissipation  rate  is  given  by 


€  =  C. 


3/4 


_  3/2 

k 


(60) 


where  the  turbulence  length  scale,  & ,  must  be  specified  consistent 
with  the  expected  turbulence  structure  in  the  two-phase  flow. 

Following  Ref.  21  the  constants  and  C  will  be  taken  as  1.0  and  0.09, 
respectively . 

In  the  present  work,  Eqs.  (55)  and  (58-60)  will  be  solved  along 
with  specified  relations  for  i  and  k 

ps 
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In  the  present  two-phase  flow  analysis  the  gas  phase  species  and 
gasified  propellant  species  mass  fractions  are  not  required.  Therefore 
in  order  to  limit  computer  requirements  the  individual  species  mass 
conservation  equations  are  not  solved,  but  rather  only  total  gas  and 
solid  phase  continuity  equations  are  solved.  Therefore,  it  is 
necessary  to  consider  transport  equations  for  the  inverse  mixture 
molecular  weight  (Z)  and  the  specific  heat  at  constant  pressure  (c  ) i 

-dy-  (apZ)  +V-(ap  UZ)  »  v[armvz]  +  Zp  r,  (61) 

where  is  the  turbulent  exchange  coefficient  for  species  diffusion 
which  is  defined  from  a  knowledge  of  the  Schmidt  number  in  the  turbu¬ 
lent  flow  of  gas  mixtures, 

r  .  -dUL  (62) 

•  m  Sceff 

and  Sc^££  is  generally  taken  as  a  constant,  Sc^^  =  0.9.  Further, 

Zp  is  the  inverse  molecular  weight  of  the  propellant  and  is  the  mass 
source  due  to  propellant  burning. 

A  similar  transport  equation  may  be  derived  for  the  specific 
heat  by  assuming  that  the  species  specific  heats  are  independent  of 
temperature : 

d(ctpC_)  r  -! 

+V-(apUcp)  =V-[armVcp]  +  (cp)pr,  (63) 


where  (0^)^  is  the  specific  heat  at  constant  pressure  of  the  propellant. 


Particle  Radius  Equation 


The  average  particle  radius,  r  ,  is  required  as  a  function  of 
spatial  location  and  time  for  the  constitutive  relations  specified 
below.  For  inviscid  flow  this  equation  may  be  written  as 
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(64) 


+  VVrp 


<  d  > 


•  2_ 

where  <d>  >0  for  surface  regression.  The  corresponding  equation 
including  turbulent  diffusion  is 


=  V.[(l-a)rmvrl-(l-a)»  (l  +  r  4L)<d> 


(65) 


where  the  relation  for  T^,  Eq.  (16),  has  been  incorporated  in  order  to 
cast  the  equation  for  the  average  particle  radius  r^  into  the  above 
form. 


Constitutive  Relations 

The  necessary  constitutive  relations  include  a  gas  phase  equation 
of  state,  a  caloric  equation  of  state,  a  turbulence  length  scale 
distribution,  the  molecular  viscosity  and  thermal  conductivity,  the 
so-called  form  functions  for  the  surface  area  and  volume  of  the  solid 
particles,  an  intergranular  stress  relation,  interphase  drag  and  heat 
transfer  relations,  and  a  burning  rate  correlation  for  the  solid  phase 
combustion.  In  the  following,  the  double  overbar  (  )  is  dropped  for 
simplicity . 

Equation  of  State  of  Gas 

The  Noble-Abel  equation  of  state  will  be  used  for  the  gas, 

p*uT 

p(l-/>7))  =  — —  =  p  ZT  (66) 

wm 

where  R  is  the  universal  gas  constant,  W  is  the  gas  molecular  weight, 
u  m 

and  q  is  the  covolume  factor,  which  is  composition  dependent. 

Following  Gough  (Ref.  1)  an  arithmetic  average  will  be  used  for  n 
based  upon  the  propellant  properties. 
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The  caloric  equation  of  state  is  taken  as 


e 


cvT 


(67) 


where  c^  is  dependent  on  the  gas  composition  but  not  temperature. 
The  static  enthalpy  is  then 


h  = 


C  T  + 
v 


ZT 


(68) 


The  specific  heat  at  constant  pressure  is 

dh  , 

cp  =  Trip  °  cv +  2 

so  that  Eq.  (68)  may  be  written  as 

h  -  cpT  +  -q p 


(69) 


(70) 


Turbulence  Length  Scale 

For  the  evaluation  phase  of  the  present  effort,  the  turbulence 
length  scale  would  be  chosen  based  upon  known  steady  state  relations. 
In  particular,  the  length  scale  would  be  taken  as  the  minimum  of  the 
length  scales  based  upon  the  local  average  distance  between  solid 
particles,  the  local  value  computed  from  turbulent  pipe  flow  correla¬ 
tions,  and  that  from  turbulent  boundary  layer  length  scale  distribu¬ 
tions  when  close  to  the  wall. 


Molecular  Viscosity,  Bulk  Viscosity,  and  Thermal  Conductivity  of  Gas 

The  molecular  viscosity  for  the  gas  is  determined  from 
Sutherland’s  law. 


where 


H- 

H-o 

110°K  for  air. 


(  T  f2  VS, 

\  T0  )  T  +  S, 


(71) 
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The  bulk  viscosity  for  the  gas  will  be  assumed  to  be  zero  at 
present, 

KB  *  0  (72) 

The  thermal  conductivity  may  be  determined  from  a  relation  similar 
to  Sutherland's  law,  e.g.. 


*  _  /  J_\3,z  To  +  Sg 

*0  '  T0  '  T  +  S2 

where  =  194°K  for  air. 

Form  Functions 


(73) 


The  surface  area  and  volume  of  particles  have  been  presented  by 
Gough  (Ref.  1)  for  a  variety  of  propellant  types.  In  the  present 
work,  spherical  propellant  grains  will  be  considered,  so 


vp  =  T  7rrp3 


sp  ■ 


(74) 


where  r^  is  the  mean  particle  radius  at  a  given  point  in  space  and 
time.  Other  propellant  types  could  easily  be  considered  within  the 
present  framework. 

Intergranular  Stress  Relation 

A  stress  relation  for  granular  propellant  has  been  given  by 
Gough  (Ref.  1),  Koo  and  Kuo  (Ref.  2),  and  Kuo,  et  al.,  (Ref.  3)  for 
the  case  when  the  average  stress  is  independent  of  the  loading 
history : 

2  ac~Q  aG 

n>  P  ( I  -  a)  a  _  c 

(75) 

0  if  a  >  ac 
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where  cx^  is  a  critical  porosity  above  which  there  is  no  direct  contact 
between  particles,  and  is  the  speed  of  sound  in  the  solid  phase 
specified  on  input.  This  relation  for  the  stress  is  obtained  by  quad¬ 
rature  from  the  following  equation  for  the  speed  of  sound,  a(a),  in  the 
solid  phase  (Refs.  1  and  2): 


d 

da 


[  ( l-a)Rp(a )] 


if 


if 


a<  a, 


a>  a. 


(76) 


Because  of  difficulties  encountered  in  obtaining  numerical  solutions 
with  implicit  representation  of  the  internal  boundaries  between 
propellant  and  gas  regions,  Gough  (Refs.  1,  17)  found  it  necessary 
to  implement  an  artificial  stress  term  by  replacing  Eq.  (76)  with 


where  is  a  "stress  attenuation  factor"  (Ref.  1)  which  must  be  speci¬ 
fied  in  an  ad  hoc  manner.  At  the  present  time  it  is  not  known  if  such 
an  approach  must  be  used  in  the  present  analysis,  however  it  could  be 
implemented  if  necessary. 

Interphase  Drag  Relation 

The  average  steady  state  interphase  drag  <F>  appearing  in  the 
momentum  equations,  Eq .  (27-28)  will  be  obtained  from  correlations 
for  nonfluidized  (packed)  regions  and  fluidized  (dispersed)  regions. 
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For  nonf luidized  regions  many  investigators  (e.g.,  Refs.  1,  6,  23) 
have  used  a  relation  deduced  from  Ergun's  (Ref.  24)  results  for  the 
pressure  drop  correlation  in  a  packed  bed  of  spheres,  i.e.. 


<  ^  >ERGUN 


pyiy 

6a2 


I50(l-a) 

ReT 


+  1.75 


(78) 


where  U_.  =  U  -  U  is  the  relative  velocity  between  the  gas  and  solid 
R  p  j  b 

particles,  and  Re^  is  the  Reynolds  number  based  on  particle  diameter 

and  relative  velocity,  i.e.. 


Re 


P 


i 

( 


H- 


(79) 


and 


Re  *  aRep  (80) 

Unfortunately,  Ergun's  correlation  is  valid  only  for  1  i  Re/(l-a)  5  4000 
and  0.4  ^  a  E  0.63,  hence  it  may  yield  erroneous  results  for  problems 
with  highly  convective  combustion  of  granular  propellants.  Recently, 

Kuo  et  al.,  (Refs.  2,  3)  have  presented  a  correlation  obtained  from 
cold-flow  resistance  measurements  under  nonf luidized ,  noncombusting 
conditions  valid  for  1  ^  Re/ (1-a)  -S  24000, 


<  F  > 


KUO 


I2arp 


Re0'87 

276.23  +  5.05 


(81) 


For  the  fluidized  region,  Koo  and  Kuo  (Ref.  2)  recommend  the 
following  correlation  obtained  from  the  expression  of  Anderssen 
(Ref.  25)  which  is  valid  for  0.003  -  Re  5  2000  and  0.45  -  a  4  1.0: 


23  . 

24  . 
25. 


Kuo,  K.K. ,  Vichnevetsky ,  R.,  ana  Summerfiela,  M. :  Theory  of  Flame 
Front  Propagation  in  Porous  Propellant  Charges  under  Confinement. 
AIM  J.  ,  Vol .  11,  No.  4,  1973,  pp.  444-451. 

Ergun,  S.:  Fluid  Flow  Through  Packed  Columns.  Chem.  Eng.  Progr., 
Vol.  48,  1952,  p.  89. 

Anderssen,  K.E.B.:  Pressure  Drop  in  Ideal  Fluidization.  Chemical 
Engineering  Science,  Vol.  15,  1961,  pp.  276-297. 
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-  i  />urIuri  f  '-a 

<f>and  =  — V5-  36zcV^r  +6cvr 


where  the  tortuosity  factor,  t^,  given  by  Ref.  2  is 


(I -a) 


for  0.45  <  a  <  O.  965 


for  0.965  <  a  <  1.0 


The  cross-section  factor  and  inertial  drag  coefficient  are 
defined  as 


c  ~  2tr{l-a)a! 


2.5  I  I- a 

t^Re/^a3-59)]1'3'-5-/ 


The  Anderssen  correlation  is  invalid  as  a  +  1,  as  noted  in  Refs.  1  and 
25,  and  a  limiting  value  must  be  imposed.  Furthermore,  the  relation 
for  C^,  Eq.  (85),  yields  -►  00  as  r^  0,  which  is  unacceptable 

1  1  n  c  r  IT  — (I  ri  oiyip  f-  -?  tts  a 

uuicod  v  v  at  lik;  l  xmc  < 

K  5 

Because  of  the  high  Reynold’s  numbers  (Re  ~  10  )  which  occur  in 

P 

typical  interior  ballistics  problems,  Gough  (Ref.  1)  has  recommended 

dropping  terms  in  the  above  correlations  containing  Re 

Gough  represented  the  interphase  drag  for  granular  propellant 

with  a  correlation  based  upon  Ergun’s  relation  for  a  settled  bed  and 

the  known  drag  coefficient  relation  for  an  isolated  sphere 
3  5 

(2x10  -  Re^  -  2x10  ) .  These  relations  were  patched  together  using 

Anderssen’ s  (Ref.  25)  correlation  between  tortuosity  and  porosity 
giving  the  following  result  (Ref.  1,  pp,  48-51): 
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1.75 


f  ■ 


1.75- 


I  -  a 


‘  a. 


0.45 


a  <ac 

ar  <  a  <  a, 

t  I 


0.3 


a  <  a  <  I 

i 


(86) 


where  is  the  settling  porosity  and  is  given  by 


a, = 


+  0.01986 


The  desired  relation  for  <;F>  is  given  by 


<  F  >  = 


PU  i 


U. 


6a' 


(87) 


/'ftjn 


Equations  (86-88)  have  been  incorporated  into  the  computer  code  in 
order  to  simplify  the  source  terms  appearing  in  the  governing  equations. 
Another  interphase  drag  correlation  could  be  incorporated  at  a  later 
date  if  warranted. 


Interphase  Heat  Transfer  Relation 

For  convective  heat  transfer  between  the  gas  and  solid  particles 
in  interior  ballistics  calculations,  numerous  correlations  have  been 
recommended  (e.g..  Refs.  1-3,  6,  23).  Gough  (Ref.  1)  advocates  the 
Gelperin-Eins tein  correlation  (Ref.  26)  for  the  interphase  heat 
transfer  with  granular  propellant  in  both  fluidized  and  nonfluidized 
regions.  The  Nusselt  number  for  this  correlation  is  (Ref.  1) 

Nup  -  2.0  +  0.4Rep'3  Pr1'3  (89) 


26.  Gelperin,  N.I.  and  Einstein,  V.G.:  Heat  Transfer  in  Fluidized 

Beds.  In  Fluidization,  edited  by  J.F.  Davidson  and  D.  Harrises, 
Academic  Press,  1971. 
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where  Pr  =  yc^/ic  and  <  is  the  gas  phase  thermal  conductivity,  Eq.  (73). 
The  heat  transfer  coefficient  in  Eq.  (49)  is  then 


h 


c 


(90) 


This  relation  is  considerably  simpler  than  the  Denton  and  Rowe-Claxton 
correlations  utilized  by  Kuo,  et  al.  (e.g..  Refs.  2,  3),  and  may  yield 
equally  reliable  predictions  in  view  of  the  large  variations  between 
experimental  data  and  the  existing  correlations  (Ref.  1). 

Finally,  the  interphase  heat  transfer  relation  required  in  the 
energy  equation  source  terra,  Eq.  (42),  is 


q  > 


•  ht(T  - 


V 


(91) 


where  h  is  given  by  Eq.  (52). 

Burning  Rate  Correlation 

The  steady  state  surface  regression  rate  (d  >  0)  is  given  by 
(e.g..  Ref.  1) 


d 


s 


=  B,  +  Bzp 


n 


(92) 


where  B^,  B9  and  n  have  known  constant  values.  The  phenomenon  of 
erosive  burning  is  assumed  to  be  an  acceleration  of  the  burning  rate 
due  to  the  influence  of  convective  heat  transfer  on  the  heat  transfer 
in  the  flame  zone.  The  Lenoir-Robillard  (Ref.  27)  regression  rate 
expression  is  utilized  for  this  effect, 


dE  =  ds  +  KEhcexp 


^E^P^E 

WJ 


(93) 


27.  Lenoir,  J.M.  and  Robillard,  G.:  A  Mathematical  Method  to  Predict 
the  Effects  of  Erosive  Burning  in  Solid-Propellant  Rockets.  Sixth 
Symposium  (International)  on  Combustion,  Combustion  Institute, 
1957,  pp.  663-667. 
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where  K^,  and  are  erosive  burning  constants,  determined  experimentally. 
The  convective  heat  transfer  coefficient  is  then  obtained  from  Eqs.  (89, 
90).  The  steady  state  burning  relation,  Eq .  (92),  has  been  incorporated 
into  the  computer  code  for  the  initial  phase  of  computations. 


Filler  Element  and  Projectile  Motion 


In  the  present  analysis  filler  elements  and  the  projectile  are 
treated  distinctly.  No  transverse  deformation  of  the  filler  elements 
is  permitted  and  elements  are  assumed  to  remain  planar;  therefore, 
a  quasi-'one-dimensional  lumped  parameter  formulation  (e.g.,  Ref.  1) 
may  be  employed  for  the  filler  elements.  The  appropriate  equations, 
which  have  been  stated  by  Gough  (Ref.  1)  are  repeated  here  for 
completeness.  It  is  assumed  that  there  are  N  filler  elements  between 
propellant  bed  and  the  base  of  the  projectile,  with  the  projectile 
denoted  as  element  (N+l) .  The  required  properties  for  each  element 
are  the  mass  (M^) ,  the  resistance  force  opposing  motion  (F^) ,  an 
internal  stress  (c^),  and  a  normal  wall  reaction  force  (Fw^)  for 
incompressible  elements  in  a  variable  area  tube.  The  cross-sectional 
area  of  each  element  is  assumed  to  be  equal  to  the  local  tube  area, 
and  the  stress  in  an  incompressible  element  is  assumed  to  be  isotropic. 

A  momentum  equation  is  then  written  for  one-half  of  element  i 
together  with  one-half  of  element  (i-1)  in  order  to  describe  the 
motion  of  the  interface  location,  z .  There  results, 


A ,  <x. 


-  *0^0- 


(94) 


—  {  M,.,  +  Mi’z,  -  A,o-i  - 


2~  (F'  +  Fm) 


-T1  W,’ 


for  2  <  i  <  N 


(95) 
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In  this  section  the  stress  is  taken  as  positive  in  tension 

following  Gough  (Ref.  1),  and  the  term  (-A  o^)  in  Eq.  (94)  is  the  force 
exerted  on  the  first  filler  element  by  the  gas  and  propellant  particles. 
The  mass  of  the  projectile  *s  assumed  to  be  corrected  for 

rotational  inertia;  if  I,  and  6  are  the  polar  moment  of  inertia, 
the  tube  diameter  and  the  angle  of  rifling,  respectively,  it  follows 


-  +  — 5-  ton20 

\  n+i  /actual  D  2 


The  normal  wall  reaction  is  given  by 


(zi+rzi 


if  element  i  is  not  incompressible 

/  dA  \  (98) 

)(t.  - ]  if  element  i  is  incompressible 

•\  dz  /, 


Constitutive  data  must  be  provided  for  the  stress  cl  for  elastic 

elements  or  for  plastic  elements  in  a  state  of  loading  (i.e., 

•  • 

z,  >  2L+^) ;  however,  for  rigid  elements  or  plastic  elements  in  a 

•  • 

state  of  unloading  (z^  £  z^+^) ,  one  has 


a*  0  ,  z.  ■  z. 


i  i+i 


Finally,  for  an  incompressible  element ,  i  ,•  one  has  the  continuity 
relation, 


(100) 
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Solution  Procedure 


The  development  of  the  MINT-G  computer  code  is  based  upon  an  axi- 
symmetric  version  of  the  highly  efficient  consistently  split,  linear¬ 
ized  block- implicit  solution  procedure  for  the  compressible  Navier- 
Stokes  equations  developed  by  Briley  and  McDonald  (Ref.  10-12),  and 
subsequently  extended  to  multi-component,  chemically  reacting, 
turbulent  flows  by  Gibeling,  McDonald  and  Briley  (Ref.  13).  This 
procedure  solves  the  Navier-Stokes  equations  written  in  primitive 
variables;  in  the  MINT-G  procedure,  the  governing  equations  are 
replaced  by  the  Crank-Nicholson  time  difference  approximation.  Terms 
involving  nonlinearities  at  the  implicit  time  level  are  linearized  by 
Taylor  series  expansion  about  the  known  time  level,  and  spatial 
difference  approximations  are  introduced.  The  result  is  a  system  of 
two-dimensional  coupled  linear  difference  equations  for  the  dependent 
variables  at  the  unknown  or  implicit  time  level.  These  equations  are 
solved  by  the  Douglas-Gunn  (Ref.  28)  procedure  for  generating  ADI 
schemes  as  per turbations  to  fundamental  implicit  difference  schemes. 
This  technique  leads  to  systems  of  one-dimensional  coupled  linear 
difference  equations  which  are  solved  by  standard  block-elimination 
methods,  with  no  iteration  required  to  compute  the  solution  for  a  given 
time  step.  An  artificial  dissipation  term  based  upon  either  a  cell 
Reynolds  number  criterion  or  the  rate  of  change  of  the  dependent 
variable  may  be  introduced  selectively  into  the  scheme  to  allow 
calculations  to  be  performed  at  high  local  values  of  the  cell  Reynolds 
number . 

The  use  of  an  implicit  solution  procedure  requires  that  equation 
coupling  and  linearization  be  considered.  Both  of  these  questions 
are  reviewed  in  detail  by  McDonald  and  Briley  (Ref.  29)  and  Briley 


28.  Douglas,  J-,  and  Gunn,  J.E.:  A  General  Formulation  of  Alternating 
Direction  Methods.  Numerische  Math.,  Vol.  6,  1964,  p.  428. 

29.  McDonald,  H. ,  and  Briley,  W.R.:  Three-Dimensional  Flow  of  a 
Viscous  or  Inviscid  Gas.  J.  Comp.  Physics,  Vol.  19,  No.  2,  1975, 
p.  150. 
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and  McDonald  (Ref.  12).  These  authors  have  argued  that  for  a  given 
grid  the  errors  arising  from  time  linearization  of  the  nonlinear  terms 
at  the  unknown  time  level  should  be  no  greater  than  the  discretization 
errors.  Also,  reduction  of  the  time  step  is  the  preferred  way  of 
reducing  the  linearization  error  since  transient  accuracy  is  thereby 
improved.  Linearization  by  Taylor  series  expansion  in  time  about 
the  known  time  level  introduces  errors  no  greater  than  those  due  to 
the  differencing  (Refs.  29  and  12),  and  this  approach  has  been  employed 
in  the  MINT-G  code.  The  formal  linearization  process  results  in  a 
system  of  coupled  equations  in  order  to  retain  second-order  temporal 
accuracy.  The  system  of  coupled  equations  at  the  implicit  time  level 
is  solved  efficiently  using  a  standard  block  elimination  matrix  inver¬ 
sion  scheme.  In  the  present  problem,  the  strong  coupling  effects 
among  the  governing  equations  dictate  the  use  of  the  block  coupled 
equation  approach.  However,  weakly  coupled  equations  would  probably 
be  solved  in  a  decoupled  manner  in  order  to  reduce  computer  time  and 
storage  requirements. 

The  principal  partial  differential  equations  which  will  be  solved 
via  the  MINT  technique  are:  gas  and  solid  phase  continuity,  gas  and 
solid  phase  momenta,  gas  phase  energy,  gas  phase  turbulence  kinetic 
energy,  gas  phase  mixture  molecular  weight  and  specific  heat  equations 
and  the  particle  radius  equation.  The  constitutive  relations  required 
to  close  the  above  system  of  equations  have  been  specified  above. 

The  solid  phase  heat  conduction  equation  is  the  only  differential 
equation  which  requires  special  treatment  because  it  is  a  Lagrangian 
equation. 

The  scheme  devised  for  solution  of  the  solid  phase  heat  conduction 
equation  is  unique  since  it  does  not  involve  the  use  of  marker  particles 
introduced  by  other  authors  (e.g..  Ref.  1).  This  is  possible  because 
the  equation  is  a  simple  heat  conduction  equation  for  a  representative 
solid  particle  moving  at  a  velocity  which  is  known  at  the  completion 
of  a  given  time  step.  The  necessary  boundary  conditions,  Eqs.  (48-52), 
provide  information  about  the  environment  through  which  the  particle 
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is  moving  in  the  form  of  a  heat  transfer  coefficient,  Eq.  (52).  The 

procedure  to  be  used  assumes  that  at  time  tn+^  the  representative 

.  ,  ,  ,  .  .  ,  -m+1  ,  n+1  n+1  n+1. 

particle  has  moved  to  the  grid  point  x.  .  ,  =  (xi .  ,  X2 .  ,  X3  )  from 

1 » J  >  K  j  k 

a  location  at  time  t  which  is  determined  from  the  known  absolute 

-►n+1  -*n  -+■ 

particle  velocities  v  and  v  ,  i.e.,  if  s  is  the  particle  position 

P  P  P 

vector  relative  to  an  inertial  reference  frame,  we  have 


as. 


dt 


(101) 


and  application  of  the  variable  time  differencing  scheme  yields, 


-n+1  -n  ,  f  .  ti  n\-n1  a* 

sp  *  sp  +  l^vP  +  (l-/3)v  J  At 


(102) 


where  8=1  for  backward  time  differencing  and  8  =  1/2  for  Crank- 

Nicholson  (centered)  time  differencing.  In  the  present  scheme, 

sn+1  is  assumed  to  be  the  grid  point  location  xn"*"^  and  Eq.  (102)  is 
P  ->n  f » J  > 

then  solved  for  s^ .  Because  the  grid  is  moving  it  is  necessary  to 

interpolate  to  find  the  required  value  of  the  particle  velocity  at 

n  .  -+n+l  -*n  -*n,->-n+l  n.  , 

time  t  at  space  point  x.  .  v  =  v  (x.  .  ,  t  ).  The  boundary 

1 » J  >  k  P  P  i  >  1  >  k 

condition,  Eq.  (51),  may  then  be  written  as 


kP  ^TP  tr  .  +  /*n*lw+n+l  +n+l\  »  ,  j./on+l  n+ul 

- - ar(C=i,t  )  =  /3[  hf(t  )  (  T  -  Tp$  )  +  kp<£(Rs  ,p  )J 

P  b 

+  (l_/3)[ht(tn){Tn  -  TpS  )  +  kp^(Rsn,pn)]  (1°3) 


The  desired  properties  h  (tn)  ,  T  n,  Tn  ,  and  cj)(Rn,pn)  are  understood 

c  p  s  s 

-fn 

to  be  evaluated  at  the  point  s^ ,  and  these  will  be  evaluated  by 
interpolation  utilizing  values  at  time  tn  at  the  four  grid  points 
surrounding  point  sn. 

Finally,  the  governing  equation  (45)  and  boundary  conditions,  Eq . 
(48)  and  (103),  may  be  written  in  finite  difference  form.  The  resulting 
tridiagonal  matrix  is  easily  inverted  using  Gaussian  elimination  to 
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yield  the  temperature  distribution  within  the  particle.  Another  approx¬ 
imate  solution  technique  could  be  incorporated  at  a  later  time  in  order 
to  reduce  the  computer  requirements  for  the  particle  heat  conduction 
model . 

Initial  and  Boundary  Conditions 


The  initial  conditions  for  the  first  phase  of  two-dimensional 


calculations  will  consist  of  a  description  of  the  fluidized  state 
of  the  flow  in  a  gun  barrel  after  ignition  is  complete  and  the  project¬ 
ile  motion  has  begun.  Typically*  the  necessary  data  would  be  produced 
from  an  existing  one-dimensional  interior  ballistics  computer  code, 
and  would  then  be  extended  over  the  two-dimensional  computational  domain 
by  applying  a  correction  for  the  wall  boundary  layers.  Provisionally, 
the  boundary  layer  integral  method  adopted  by  Gough  (Ref ,  1)  would 
be  utilized  to  determine  the  boundary  layer  thickness  and  velocity 
prof ile. 

The  boundary  conditions  to  be  applied  would  be  no-slip  wall 
velocities  on  solid  surfaces  and  conventional  symmetry  conditions  at 
the  tube  centerline.  The  breech  would  be  assumed  to  be  stationary, 
and,  of  course,  the  projectile  and  filler  elements  would  be  allowed 
to  move.  The  wall  pressure  would  be  determined  by  employing  the 
normal  gas  momentum  equation  written  at  the  wall.  The  surface  tempera¬ 
ture  would  be  determined  by  incorporating  a  barrel  heat  conduction 
model  coupled  to  the  gas  heat  transfer  at  the  wall.  For  simplicity, 

u  4-  ^  ^  ~  *-u~  "u  ~  ~  i  - 1  A  - - * A  +-  ~  i ~  ~  i  , .  -*  ~ 

ucaL  GuuuQGLiuu  _Lu  Lilt:  udiiti  wuuiu  ue  dbbuu  leu  lu  ue  pjL  xiua  jl  x  ±y  xn 

the  radial  direction.  The  porosity  at  a  wall  would  be  determined  from 
either  the  solid  phase  continuity  equation,  Eq .  (8),  or  the  solid 
phase  momentum  equation,  Eq.  (28),  written  at  the  wall. 

The  appropriate  boundary  condition  for  the  inverse  mixture 

molecular  weight  (Z) ,  specific  heat  (c  ) ,  and  particle  radius  (r  ) 

P  P 

at  a  solid  wall  is  zero  normal  derivative,  i.e., 


(104) 
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This  follows  from  the  definition  of  these  quantities  as  mass  weighted 
averages,  and  the  assumption  that  the  individual  species  diffusion 
velocities  normal  to  the  wall  as  determined  from  Fick's  law  must  be 
zero;  that  is,  (9m^/3n)  =  0  where  itu  is  a  species  mass  fraction. 

The  solid  particles  which  reach  the  wall  will  be  assumed  to  be 
in  equilibrium  with  the  gas  phase,  thus 

(  =  (T)„  (105) 


40 


THE  COORDINATE  SYSTEM 


The  set  of  governing  partial  differential  equations  which  model 
the  physical  processes  occurring  in  interior  ballistics  problems  was 
presented  in  the  previous  section.  For  generality  these  equations 
were  written  in  vector  notation;  however,  before  these  equations  can 
be  incorporated  into  a  computer  code,  a  coordinate  system  must  be 
chosen.  The  governing  equations  can  then  be  cast  in  a  form  reflecting 
the  choice  of  the  coordinate  system.  In  choosing  a  coordinate  system 
for  interior  ballistics  calculations,  it  was  felt  that  there  are  two 
primary  considerations:  (1)  the  coordinate  system  must  have  the  ability 
to  enlarge  the  physical  extent  of  the  computational  domain  as  the  pro¬ 
jectile  moves  through  the  gun  barrel,  and  (2)  the  coordinate  system 
must  be  of  a  general  enough  nature  such  that  future  modifications  to  the 
geometry  portion  of  the  computer  code  can  be  accomplished  without  a 
major  restructuring  of  the  code.  With  the  above  in  mind,  it  was 
decided  to  utilize  a  moving  three-dimensional  general  orthogonal 
coordinate  system.  The  governing  equations  were  obtained  by  special¬ 
izing  the  moving  three-dimensional  general  nonorthogonal  equations 
presented  by  Walkden  (Ref.  30)  for  the  present  interior  ballistics 
problem.  Because  the  equations  of  Walkden  consider  a  moving  coordinate 
system  consideration  (1)  above  is  satisfied.  In  addition  it  is  felt 
that  consideration  (2)  is  satisfied  by  the  use  of  three-dimensional 
general  orthogonal  coordinates,  since  the  geometries  associates  with 
most  interior  ballistic  problems  can  adequately  be  described  by  such 
a  system.  The  set  of  governing  partial  differential  equations  for  this 
coordinate  system  is  presented  in  Appendix  A.  In  these  equations  h^, 
h2  and  h^  represent  the  metric  coefficients  in  the  x^,  and  x^ 
coordinate  directions,  respectively,  and  the  Jacobian,  J,  is  defined 

by 

J  =  h,h2h3  (106) 

30.  Walkden,  F. :  The  Equations  of  Motion  of  a  Viscous,  Compressible 
Gas  Referred  to  an  Arbitrarily  Moving  Coordinate  System.  Royal 
Aircraft  Establishment,  Technical  Report  No.  66140,  April  1966. 
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The  gaseous  velocity  components  are  represented  by  u,  v  and  w  for  the 
Xl*  X2  an<^  X3  coordinate  directions,  respectively,  while  the  corres¬ 
ponding  solid  phase  velocity  components  are  u^,  v^  and  w^.  It  is 
assumed  that  the  projectile  moves  only  in  the  x^-coordinate  direction 
and  hence  the  metric  has  the  functional  form 


h3  *  h3(X3«,) 

In  addition  it  is  assumed  that 


(107) 


and  that 


h|  -  constant 
h2  =  h2(x() 


(108) 


(109) 


For  example,  in  cylindrical  polar  coordinates,  h^  -  1  and  h^  =  x^.  The 
last  two  assumptions  considerably  simplify  the  analysis  of  the  viscous 
stress  terms  and  apply  strictly  for  the  coordinate  systems  of  interest 
under  this  effort. 

In  the  present  formulation  the  x^-direction  velocity  components 
(w  and  Wp)  are  measured  relative  to  the  moving  coordinate  system.  The 
terms  which  include  the  effect  of  the  moving  coordinate  system  appear 
only  in  the  x^-direction  momentum  equation  (since  the  projectile  motion 
is  limited  to  the  x^-direction) .  In  the  x^-direction  gas  phase  momentum 
equation  this  results  in  the  additional  time  term 


.  .  ^Vg 

h3 JaP  dt 

while  the  viscous  stress  term  is  augmented  by 


(110) 


2h. 


Re  dx 


a  M 


h2  dvQ 


dx. 


(Ill) 


where  v  is  the  grid  velocity  (in  the  x„-direc tion) . 

By  definition  the  x^-direc tion  metric  can  be  expressed  as 


h,(x_  ,t)  = 


dz. 


dx 


(112) 
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where  z is  the  cartesian  (physical)  coordinate.  In  order  to  allow 
for  a  nonuniform  physical  grid  in  the  x^-direction,  a  transformed 
normalized  coordinate,  n(x^)  is  defined  by 


t7(x3) 


(113) 


where  zq  is  the  cartesian  location  of  the  breech  end  of  the  gun  barrel 
and  z i  is  the  cartesian  location  of  the  first  filler  element  end  and 
x^  is  now  specified  as  being  equally  spaced  and  having  a  value  from 
1.0  to  Combining  Eqs.  (Ill)  and  (112)  yields 


h_(x_,t)  =  (z-  z J 


drj 


o'  dx_ 

3 


(114) 


Note  that  the  time  dependence  of  h^  is  introduced  through  z^  which 
varies  as  the  projectile  moves  through  the  gun  barrel.  The  local  grid 
velocity,  v  ,  can  be  calculated  from  the  relationship 


(115) 


where  z^  represents  the  velocity  of  the  first  filler  element. 

The  functional  form  of  n(x^)  is  arbitrary  and  can  be  chosen  such 
that  the  packing  of  grid  points  in  the  x^-direction  is  achieved  in  the 
regions  where  the  largest  gradients  are  expected.  Presently  the  compu¬ 
ter  code  allows  for  the  concentration  of  gr-id  points  to  occur  by  means 
of  a  generalization  of  the  Roberts’  transformation  (Ref.  31).  The  grid 
points  can  be  concentrated  at  the  breech  end  of  the  computational  domain 
and/or  at  the  filler  element  end  of  the  computational  domain  or  the  grid 
points  can  be  concentrated  around  some  interior  location.  The  trans¬ 
formation  equation  used  for  this  purpose  is 


V 


/ 


C 


tanh(Ex3  +  F)-H 
_  -  — 


+  D 


(116) 


31.  Roberts,  G.E.:  Computational  Meshes  for  Boundary  Layer  Problems. 
Proceedings  of  the  Second  International  Conference  on  Numerical 
Methods  in  Fluid  Dynamics,  Springer-Verlag ,  New  York,  1971,  p.  171. 
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where  no  is  the  value  of  n  about  which  the  concentration  of  grid  points 
is  centered  and  the  values  of  A,  C,  D,  E,  E,  G  and  H  are  controlled  by 
the  input  parameters,  r)Q,  t^,  and  The  derivation  of  the 

relationships  between  A,  C,  D,  E,  F,  G  and  H  and  the  input  parameters 
is  quite  lengthy  and  hence  only  the  results  are  presented  here,  viz.. 


sinh(t2) 
’W  '  Vq 


(in  this  study  =  1.0) 

MAX 


where 


C  - 


m  7^ 


(118) 


t,  5  ln{A<VMIN  -Vo )  +  V~'  +  a2  (  VM(N- Vo)2  }  (119) 


(in  this  study  nMIN  =  0) 


D  =  t,  -  C 


(120) 


where 


(121) 


(122) 


(123) 


(124) 


44 


G 


(125) 


-  I 


MAX 


H  =  T.  -  G  (126) 

The  above  is  presented  only  for  completeness;  the  important  thing  to 
note  is  the  effect  that  r^,  and  t ^  have  on  the  physical  grid 

SDacine.  The  effect  of  t~  is  to  reeulate  the  sinh  oortion  of  the 
transformation,  while  x^  and  regulate  the  tanh  portion  of  the 

transformation;  x^  controls  the  physical  grid  spacing  at  the  breech 
end  of  the  computational  domain  while  x ^  controls  the  spacing  at  the 
first  filler  element  end.  The  values  of  x^  and  x^  are  subject  to  the 
following  limitations 


-  1  <  T,  <  0 

(127) 

0  <  Tg  <  1 

(128) 

ro 

tv 

O 

(129) 

In  order  to  see  how  the  input  parameters  effect  the  grid  spacing  it  is 

instructive  to  first  negate  the  effect  of  the  sinh  by  setting  =  0 

and  investigating  the  effect  that  x^  and  x^  have  on  the  transformation. 

If  x^  =  0  and  x^  >  0  grid  packing  will  occur  at  (the  larger  the 

value  of  x 2  the  greater  the  packing)  while  if  x^  <  0  and  x^  =  0 

packing  occurs  at  (the  larger  the  value  of  |x^|  the  greater,  the 

packing) .  Zero  values  of  x^  and  x^  result  in  equal  grid  spacing  while 

nonzero  values  of  both  x^  and  x ^  result  in  packing  at  both  p  ^  ana 

^MAX*  °^-^er  han<3  11  the  effect  of  the  tanh  is  negated  (by 

setting  both  x^  and  x^  both  equal  to  zero)  the  effect  is  to  concentrate 

the  grid  points  about  n  only.  The  larger  the  value  of  t^  the  greater 
-  -  O'  _  z 

the  concentration.  Nonzero  values  of  t2>  x^  and  x^  result  in  a 
combination  of  the  effects  of  the  sinh  and  the  tanh  transformations. 
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Equation  (116)  allows  one  to  concentrate  the  physical  location  of 
grid  points  in  the  x^-direction  in  some  prescribed  manner.  However 
it  is  also  desirable  to  have  a  transformation  technique  which  concen¬ 
trates  grid  points  in  the  x^-direction  as  a  function  of  and  time,  t. 
Such  a  technique,  for  instance,  would  permit  the  concentration  of 
x^-direction  grid  points  in  the  manner  shown  in  Fig.  2  to  account  for 
a  variation  of  the  boundary  layer  thickness  as  a  function  of  x^  and  t 
as  the  projectile  moves  through  the  barrel.  As  can  be  seen  in  Fig.  2 
the  resulting  coordinate  system  is  nonorthogonal ,  but  since  the 
nonorthogonality  is,  so  to  speak,  only  in  the  x^-direction,  the 
increased  degree  of  computational  difficulty  is  not  large.  The 
partial  differential  equations  in  coordinate  systems  like  that  shown  in 
Fig.  2  can  be  obtained  by  transforming  the  orthogonal  governing 
equations  in  Appendix  A.  In  order  to  allow  the  greatest  degree  of 
flexibility  only  a  general  functional  form  of  the  transformed  variable, 
y^,  will  be  prescribed  at  this  time,  viz., 

y ,  ■  y,(x|fx3I t)  030) 

Therefore  a  general  variable,  <j) ,  is  transformed  by  the  relationship 

<£(x,,x2lx3lt)  -  <#>(y|Iyg)y3,t)  U3i) 


where 

(132) 

(133) 

By  use  of  the  chain  rule  the  space  and  time  derivatives  of  <{>  can  be 
calculated.  For  example  the  calculation  of  9<fj/3xj_  proceeds  as  follows: 


y2  =  x2 

y3  =  s 


d<f>  dcp  dy 


dx.  dy(  dx (  ’  dy2  '  dy3  .  '  dt"  (13A) 


dc£  dt/ 


0  =0  =0 

;rni 

are  replaced  by  first  derivatives  with  respect  to  y^  times  a  scaling 
factor  3y^/9x^  (which  can  be  calculated  once  the  functional  form  of 
Eq.  (130)  is  prescribed).  The  calculation  with  respect  to  x„ 


Thus  the  first  derivatives  with  respect  to  x^  in  the  governing  equations 


proceeds  as  follows: 
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defy 


dx 


3 


defy  dy  j 


(135) 


The  first  derivative  in  the  x^-direction  can  thus  be  replaced  by  the 
first  derivative  in  the  y^-direction  plus  the  first  derivative  in  the 
y^-direction  times  the  scaling  factor  8y^/9x^.  It  can  easily  be  shown 
that  first  (and  second)  derivatives  in  the  X2~direction  have  the  same 
form,  i.e., 


defy  d<p 

J*l  ~  ~ay \ 


(136) 


The  time  derivative  of  <£  is  calculated  from 


d<f> 

JT 


defy  dy, 

d^  dT 


(137) 


Therefore  the  time  derivatives  are  replaced  by  a  time  derivative  plus 
a  convective  like  term  times  a  scale  factor  9y^/8t.  Further  use  of  the 
chain  rule  yields  the  equations  for  second  derivatives.  They  are 


dzef>  d2y,  d<fy  /  dy,  \2  dzef> 
dx,4  *  dx,2  dy7  +  \’ax7/  ay, 2 

dz<£  dz<f> 

ax22  =  ay22 


(138) 


(139) 


d*4>  dz4>  /dy,  \2  d2<£  a2 y ,  dcf> 
dx33  =  dy32  +  \  dx3  )  ay,2  +  dx32  dy. 


dy,  d*<f> 

+  Z~d^~d^dTl 


(140) 


and 


d2<£ 


dx,dx3 


dy, 

dx. 


dz<p  dy. 


d2</> 


dy,2  dx3 


ayt  aV3. 


defy 

dy,  dy,  dy. 


azy, 


(i4i) 
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Equations  (134)  through  (141)  are  used  to  transform  the  orthogonal 
version  of  the  governing  equations  into  their  nonorthogonal  form.  The 
results  are  presented  in  Appendix  B.  It  is  important  to  note  that 
although  the  derivative  operators  have  been  transformed  into  a  non¬ 
orthogonal  frame,  the  velocity  components  are  still  the  original 
orthogonal  (x^  x2>  x^)  components  of  velocity  with  the  x^direction 
velocity  components  (w  and  w^)  measured  relative  to  the  coordinate 

system  motion  at  the  grid  velocity  (v  )  in  the  x„-direc tion . 

g  3 

The  specific  functional  form  of  Eq .  (130)  has  not  yet  been 
specified.  From  the  governing  equations  of  Appendix  B  it  can  be  seen 
that  the  following  scale  factors  are  needed,  viz.,  3y../3x  ,  3y../3x~, 

2  2  2  2  n  -L-L-LJ 

3  y^/3x^  ,  3  y^/3x^  ,  3  anc*  3y^/3t.  The  computer  code  is 

set  up  in  such  a  manner  that  the  user  can  implement  a  variety  of 
functional  forms,  hence  only  a  specific  functional  form  is  presented 
here.  The  form  chosen  here  is  a  generalized  Roberts’  transformation 
with  coefficients  that  vary  with  x^  and  t,  i.e., 

Ey (  +  F  r  tanh"1  { Gx(  +  H )  (142) 

In  the  above  equation  the  parameters  E,  F,  G  and  H  are  defined  in  a 
similar  manner  as  in  Eq .  (116).  The  relationships  are 

tanh" 't  - tanhHr. 

E  =  - - - ; - 1  (143) 

y.  - 1 

MAX 


y  tanhf'r  -  tanh_lT 

max _ [ _ 

y.  - 1 

MAX 


(144) 


I 

MAX  'MIN 


(145) 
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(146) 


H  = 


x 


maxT|  X|minTz 

xl  _xl 
MAX  MIN 


where  for  the  interior  ballistics  problem  refers  to  the  axis 

of  symmetry  and  thus  =  0*  Thus  as  before  these  parameters  are 

controlled  by  values  of  and  for  the  x^-direction.  The  difference 

is  that  in  this  case  the  and  T2  can  be  functions  of  x^  and  t. 

Eq .  (142)  has  the  capability  of  forming  an  adaptive  mesh  generator 
which,  for  example,  could  follow  the  sidewall  gun  barrel  boundary 
layer  growth  as  it  develops  as  a  function  of  both  and  t.  By 

differentiating  Eq.  (142)  with  respect  to  the  proper  variables  all  of 
the  above  scale  factors  can  be  calculated.  For  example  differentiation 
of  Eq.  (142)  with  respect  to  x^  yields 


dy,  g  1 

ax,  =  i-tz 


(147) 


where 


T  s6x|  +  H 


(148) 


and  likewise 

a2y,  2GZT 

ax(2  '  e(i-t2)2 


(149) 


Differentiation  with  respect  to  x^  and  t  becomes  slightly  more  compli¬ 
cated  as  and  (and  thus  E,  F,  G  and  H)  are  functions  of  and  t. 
For  example  differentiation  with  respect  to  t  yields 
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Differentiation  with  respect  to  x~  yields  the  same  formula  as  above, 

J  2  2 

except  that  t  is  replaced  with  x^.  The  formula  for  9  y^/^x^  is  quite 
lengthy  and  as  might  be  expected,  involves  both  first  and  second 
derivatives  of  the  and  with  respect  to  x^,  i.e., 


dr 


a  y, 

ax32 


aX32 


dr. 


r  U, 


dr. 


2T 


"V 

MAX  3 


+U|'X'MlN'aX3 


L 

,  "xi 

MAX 


(x,  —  x,  )(l-T2)2 

IN 
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‘max  'min 


MAX 


-  X 


MAX 


r  xrx, 


_  x  2 
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MIN 


MAX 


dx. 


MIN 


X  -  X 
L  MAX  MIN 


y, -i 


Ft* 


-  T 


MAX 


y,  -i 

MAX 


I  dr, 


l“T,2  *X3  L 


dy 


dx 


+  ( y,  -y,)r 


2  'J 

dr 


MAX 


i'  » ax 


3  J 


I  dr. 


I  - T  2  dX 
2  *'"3 


ay,  a-r. 

— +  (y.-')T 


dx 

-"3 


3  J 


(151) 
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For  the  functional  form  given  by  Eq -  (142),  the  formulae  for  the 
scaling  factors  require  a  priori  knowledge  of  their  first 

derivatives  with  respect  to  x^  and  t  and  their  second  derivatives 
with  respect  to  x^  in  order  to  make  the  mesh  adaptive.  Values  of 

and  ^2  an<^  their  space  and  time  derivatives  can  be  specified  by  two 
techniques.  In  the  first  technique  the  functional  form  of  i^(x^,t) 
and  i2(x2,t)  would  be  specified  and  then  differentiated.  Preferably 
the  functional  form  would  be  of  such  a  nature  so  as  to  concentrate 
grid  points  in  the  areas  of  largest  gradients.  The  second  technique 
would  involve  scanning  the  solution  field  at  time  t  for  the  areas 
of  large  gradients,  and  choosing  T^Cx^t)  and  T2(x^,t)  such  that 
grid  point  packing  occurred  in  those  areas.  Values  of  T^(x^,t)  and 
T2^X3,t:^  could  then  be  stored  and  their  derivatives  calculated 
numerically. 

In  summary,  this  section  has  presented  the  techniques  for  taking 
the  vector  form  of  the  governing  partial  differential  equations  which 
model  the  two-phase  flow  phenomena  occurring  in  interior  ballistics 
problems  and  has  converted  them  into  a  form  which  can  be  directly 
programmed  into  a  computer  code.  The  resulting  equations  take  into 

consideration  the  expansion  of  the  computational  domain  as  the  project¬ 
ile  moves  through  the  gun  barrel.  Rather  general  transformation 
techniques  are  used  to  concentrate  grid  points  in  areas  of  steep 
gradients  in  both  the  axial  and  radial  directions.  The  resulting 
equations  (see  Appendix  B)  are  nonorthogonal ;  however,  the  velocity 
components  are  still  in  the  original  orthogonal  (x^,X2»x^)  coordinate 
directions  with  the  x^-direction  velocity  components  (w  and  w  ) 
measured  relative  to  the  coordinate  system  motion  at  the  grid  velocity 
(v  )  in  the  x--direc tion . 

s  ^ 
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APPENDIX  A 
GOVERNING  EQUATIONS 


Gas  Phase  Continuity  Equation 


—  (Ja/S)  .  -{  ±-  (  h2h3a^U) 

d  _  a 


+  <h.h3a^ v)  +  ir(hih2a^w)} +  jri 


(A-l) 


Solid  Phase  Continuity  Equation 


+  rtW^Vp]  +  ^-[h^d-a^pW  ]}  -  jr, 


(A-2) 


Gas  Phase  Momentum  (x^-Direc tion) 
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+  h .  J 


[  -  ( I  -  a)  — —  <F>!  +  u  r 

L  v  I  P  i. 


59 


Solid  Phase  Momentum  (x^-Direc tion) 
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Gas  Phase  Momentum  (x2-Direction) 
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Solid  Phase  Momentum  (^-Direction) 
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Gas  Phase  Momentum  (x3~Direction) 
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Solid  Phase  Momentum  (x3-Direction) 


at 
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Gas  Phase  Energy  Equation 
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and  A  ^  is  the  nondimensional  form  of  Eq.  (42), 


a,  *( — — )(e  - 0- m  +  r  ) 
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General  Equation 
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where  tp  represents  either  the  inverse  of  gaseous  mixture  molecular 
weight  or  the  gaseous  mixture  specific  heat  and  <j>p  represents  the 
corresponding  propellant  property. 

Particle  Radius  Equation 
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APPENDIX  B 

GOVERNING  EQUATIONS  (TRANSFORMED) 
Gas  Phase  Continuity  Equation 
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Solid  Phase  Continuity  'Equation 
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(B-6) 
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Solid  Phase  Momentum  (xj-Direc tion) 
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Gas  Phase  Momentum  (x2-Direction) 
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Gas  Phase  Momentum  (x3~Direction) 
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Solid  Phase  Momentum  (x-j-Direc  tion) 
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Gas  Phase  Energy  Equation 
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where 
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and  A  ^  is  given  in  Eq .  (A-ll) . 


General  Equation 

d 

—  {Jap<j>)  *  - 
dy,  d 


dy  d  dy.  d 

777  777  +  777  777  ‘  h,h2a^> 


d  d 

dT  ~UaP^]  +  77"  (  h,h3a^v^)  +  —  (h,h2a^w<#>) 


(B-23) 


ReSc 


q  d26 

G_  +  G-  _  . 

"3  rj  dy« 


7T  +  [g.+G.+G_+G-1  "T^ —  + 


L  I 


d<£ 

5  ^bj  ay 


h  I  h3  d 


/  d4>  \ 

I  au: - I 


h0  0yo  \  '  dy,/ 


h,h2 


d24> 


+  G 


d<£  2h,h2  ay,  d2<£ 


2  ay. 


dx,  ^  dy,dy. 


+  ^Pri 


72 


where  4>  represents  either  the  inverse  of  gaseous  mixture  molecular 
weight  or  the  gaseous  mixture  specific  heat  and  <J>p  represents  the 
corresponding  propellant  property. 

Particle  Radius  Equation 
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